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Abstract 

We examine the effect of the phase-resetting curve (PRC) on the 
transfer of correlated input signals into correlated output spikes in 
a class of neural models receiving noisy, super-threshold stimulation. 
We use linear response theory to approximate the spike correlation 
coefficient in terms of moments of the associated exit time problem, 
and contrast the results for Type I vs. Type II models and across the 
different timescales over which spike correlations can be assessed. We 
find that, on long timescales, Type I oscillators transfer correlations 
much more efficiently than Type II oscillators. On short timescales this 
trend reverses, with the relative efficiency switching at a timescale that 
depends on the mean and standard deviation of input currents. This 
switch occurs over timescales that could be exploited by downstream 
circuits. 



1 Introduction 

Throughout the nervous system, neurons produce spike trains that are cor- 
related from cell-to-cell. This correlation, or synchrony, has received ma- 
jor interest because of its impact on how neural populations encode in- 
formation. For example, correlations can strongly limit the fidelity of a 
neural code as measured by the signal-to-noise ratio of homogeneous pop- 
ulations [HI EH El II]- However, the presence or stimulus-dependence of 
correlations can also enhance coding strategies that rely on discriminating 
among competing populations [21EJES]; in general, the effects of correlation 
on coding are complex and can be surprisingly strong [33 U2] EE E3 13 
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SH [261 El 113 HZ1 SHI S3]. In addition, stimulus-dependent correlations can 
modulate or directly carry information directly |15| S3 ESI El EU El [HI El] • 
Correlations also play a major role in how signals are transmitted from 
layer-to-layer in the brain [JTJ E3 [28] . 

What is the origin of correlated spiking? One natural mechanism is the 
overlap in the inputs to different neurons - these common inputs can drive 
common output spikes. This poses the question: how does the process of 
transferring of input correlations to spike train correlations depend on the 
nonlinear dynamics of individual neurons? Such correlation transfer has 
been modeled in integrate-and-fire type neurons jH [32l Q31 SHJ HQ [28j and, 
very recently, in phase reductions of neural oscillators |31| ITT?] . 

In particular, |31[ Qjjj contrast the correlated activity evoked in neural 
oscillators with Type I (i.e., always positive) vs. Type II (i.e., positive and 
negative) PRCs. When correlations are measured via equilibrium probabil- 
ity distributions of pairs of neuron phases, or via cross-correlation functions 
of these phases over time, Type II oscillators display relatively higher levels 
of correlation |31[ [T9] . These results for phase correlation imply a similar 
finding for spike train correlation in the limit of very short timescales (the 
connection arises because the spike train cross-correlation function at r = 
can be related to the probability that phases will be nearly coincident for 
the two cells {35].) 

In this study, we also contrast correlation transfer in Type I and Type II 
oscillators (as well as in a continuum of intermediate models). The primary 
extension that we make is to study spike train correlation over a range of 
different timescales. Specifically, we measure the correlation coefficient px 
between the number of spikes produced by a pair of neurons in a time window 
of length T: 

p T - Cov(n 1 ,n 2 ) _ ^ 

y / Var(n\)y / Var(n2) 

Here, m, ri2 are the numbers of spikes output by neurons 1 and 2 respectively 
in the time window; see Fig. 1 for an illustration. 

We first derive a tractable expression for px in the long timescale limit 
T — > oo (cf . |144 SB] ) . This can be given in terms of moments of an associated 
exit time problem. This reveals a dramatically higher level of long-timescale 
correlation transfer in Type I vs. Type II neural oscillator models, the 
opposite of what was found in the earlier studies over short timescales (see 
Fig. 1). Next, we study pt for successively shorter T, recovering the earlier 
findings of |19| [3~T] . and noting the critical timescale below which Type 
II neurons become more efficient at transferring correlations. Additional 
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Figure 1: (a) A schematic of the setup and correlation metric used in this 
study, (b) The principal result of our study; that correlation transfer effi- 
ciency depends on both internal dynamics and timescale of readout, with 
relative efficiency between Type I (light red) and Type II (dark blue) switch- 
ing as readout timescale changes. The graph shows pT for a particular set 
of model parameters vs. logarithm of the time window log(T). Insets show 
the phase-resetting curves used for the Type I and Type II models. 




results on how correlation transfer depends on neurons' operating range - 
that is, their spike rate and coefficient of variation (CV) of spiking - are 
developed as we go along. 

2 Models of neural oscillators and correlation trans- 
fer 

2.1 Phase reductions 

Neural oscillators can be classified into two types based on their intrinsic 
dynamics. Both types have the feature that as applied inputs (or "injected 
current") increases, the system transitions from a stable rest state to periodic 
firing through a bifurcation, the nature of which defines the type [38]. Here, 
as is often taken to be the case, Type I neurons undergo a saddle-node on 
invariant circle bifurcation, in which two fixed points collide and disappear, 
producing a periodic orbit that can have arbitrarily low frequency. Type 
II neurons undergo either a subcritical or supercritical Hopf bifurcation, in 
which a periodic orbit emerges at a non-zero minimum frequency. 

Once the oscillator has passed through this bifurcation, it can be de- 
scribed by a single equation for its phase, in which inputs are mediated 
through a phase-resetting curve (PRC) which indicates the degree to which 
an input advances or delays the next spike. A PRC is derived for a mathe- 
matical model by phase reduction methods and determined experimentally 
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by repeated perturbations of a system by input "kicks" [50 . Several in- 
vestigators have demonstrated a connection between the type of bifurcation 
and the shape of the PRC: a Type I oscillator PRC is everywhere positive, 
so that positive injected current advances the time of the next spike |17j . 
whereas a Type II PRC has both positive and negative regions, so positive 
inputs advance or delay the next spike depending on their timing ( |18j : see 
also [9]). Moreover, the form of Type I and Type II phase-resetting curves 
near the bifurcation are given by (1 — cos0) and — sin(0) respectively. We 
investigate these two PRCs, together with a family of PRCs given by a linear 
combination of these prototypical examples: 

Z(0) = -asin(0) + (1 -a)(l -cos(0)), < a < 1 . (2) 

Here a homotopes the PRC from "purely" Type I to Type II; along the way, 
intermediate PRCs more representative of phase reductions of biophysical 
models are encountered [T71 [22] . 

The question of how oscillators with different PRCs synchronize when 
they are coupled has been the subject of extensive study. Here, we ask about 
a different mechanisms by which such oscillators can become correlated. 
Specifically, we consider an uncoupled pair of neurons receiving partially 
correlated noise. The dynamics have been reduced to a phase oscillator; 
each neuron is represented by a phase only. Each phase 9i, is governed by 
the stochastic differential equation 

al9 i = ujdt + (jZ{e i )o{^l^~cdWi + ^cdW^), 0iG(O,2vr) (3) 

where uj, a > 0, o denotes the Stratonovich integral, and 

Z{6) = -asin(0) + (1 - a)(l - cos(0)) . 

Each 6i receives independent white noise, dW\, and common white noise 
dWf is received by both. The noises are weighted so that the total variance 
of the noise terms in (J3j) is a 2 . For the remainder of the paper, we treat the 
equivalent ltd integral 

2 

d6i = (uj + ^-Z{0i)Z'{9i)) dt + aZ(9i)dW t , 9 { £ (0, 2tt). (4) 
2.2 Measuring spike train correlation 

We record spike times as those times t\ at which 9\ crosses 2tt [T71 [TS] . 
Because Z{2ir) = and uj > for all the models we consider, 9 always 
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continues through 2ir and begins the next period of inter-spike dynamics. 
We consider the output spike trains yi(t) = ^ S(t — i*), where t\ is the time 
of the fcth spike of the ith neuron. The firing rate of the ith cell, (yi{t)), is 
denoted V{. As a quantitative measure of correlation over a given time scale 
T, we compute the following statistic: 

Cov(ni,n 2 ) 

PT ~ 



yVar(n\) \JVar(n 2 ) 

where n\ , n 2 are the numbers of spikes output by neurons 1 and 2 respec- 
tively in a time window of length T; i.e. rii(t) = // +T Vi(s) ds. 
One can show that this is equivalent to 

f T T C 12 (t)^dt 
PT = , (5) 
Xt (t) ^ dt /5 T C 22 (t) ^ dt 

where Cij(r) = (yi(t)yj(t + r)) — ViVj |13j . It will be convenient for us to 
analyze the system in the Fourier domain. By the Wiener- Khinchin theorem, 
we can write ^ in terms of the power spectra Pij = (y*yj) as 

IZoPMK T (f)df 
PT = . (6) 
y/fT^ Pn(f)K T (f) df P 22 (f)K T (f) df 



where the kernel Kt is 



KM) = 



3 Correlation in the long timescale limit 
3.1 Linear Response Theory for px 

We recall the following derivation from [29] HH H5] . Assume that the fraction 
of noise variance c of the correlated input noise is small; then we treat the 
system with common noise as a perturbation to the system without common 
noise. Because the common noise is small, we will assume that the response 
of the system can be treated by linear response theory; that is in the Fourier 
domain it can be characterized by a susceptibility function A u a(f) which 
gives the scaling factor between input and response at frequency /. 
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Specifically, we make the ansatz [29 that the Fourier transform can be 
written to lowest order in c 



Uf) = M/) + v^W(/)Q(/) 

where yo,i is the spike output of the neuron without correlated noise, Q(f) 
is the Fourier transform of the correlated noise, and A is the susceptibility 
function. Then the cross-spectrum of the two spike trains, Pi2(/) ; satisfies 

Pn(f) = c\A^ a (f)\ 2 (Q*Q) 
= ca 2 \A^(f)\ 2 

as the base spike trains are independent of each other and the correlated 
noise Q, and taking the variance of Q to be a 2 . Then at any finite T, px is 
linear in c, and (cf. [SH], IS]): 

Pt(w,ct) = cS T (u;,cr) (7) 

^ 2 IZ\^Af)\ 2 KT(f)df (8) 
JZ Pn(f)K T (f)dfJZ P22(f)K T (f)df ' 

Note that St(u, a) multiplies the input correlation c to yield the spike train 
correlation; for this reason we refer to St(u,(j) as the correlation gain. We 
recover a simple expression for ^ in the limit T — > oo. In the numerator, 
we converge to the value of the integrand at as / — > oo: 

ca 2 \A^ a (0)\ 2 

As Aw ; a{0) is the limit of the susceptibility function as the frequency becomes 
arbitrarily small, it must be equivalent to the ratio of the DC response of 
the system (that is, the firing rate v) to the strength of a constant DC 
input. Later we will use the symbol p to include a DC input explicitly for 
the purposes of this computation. Therefore we define 

The denominator converges to -Pn(O) (assuming the unperturbed oscillators 
to be statistically identical) which for renewal processes is simply CV 2 v. 

Putting these results together, as T — > oo, the finite-time correlations 
satisfy 

lirn^ p T = c-^^y = cS(uj, a) (10) 
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where v is the mean output firing rate, CV the coefficient of variation of 
the interspike intervals, a the input noise amplitude. Each quantity can be 
computed from statistics of a single oscillator, and combined to yield the 
long-time correlation gain S[u,o). 



3.2 Computing moments of the exit time problem 

The quantities CV, and ^ are related by moments of the interspike 
intervals (ISI), and as such can be computed using the associated exit time 
problem. Specifically, 

CV = V^-f^ (12) 
Ti(0) V ; 

dv 1 dT^O) 

dn (Ti(0)) 2 dn 1 ' 

where Ti(0) is the average time to cross 8 = 2ir, given a start at 9q = (in 
other words, given a start at the last spike), and T2(0) is the second moment 
of this same quantity. 

T\(x) and T2(x) are given by solutions of the adjoint equation [50] 

fi(x) = A(x)d x T i + ^B(x)d 2 x T i (14) 

where 

fi(x) = -1 (15) 

f 2 (x) = -2Tx(x) (16) 

2 

A{x) = cj + ^-Z(x)Z\x) (17) 
B(x) = a 2 Z(x) 2 (18) 

and boundary conditions are given by Tj(2-7r) = 0, ^ bounded at x = 0. 
The solution can be obtained by integrating equation (14) twice over the 
range [0,27r], using the integrating factor ^f: 

= ^(/^^y)' *(0) = 0. (19) 
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If a 



0, then B{x) has no interior zero and we proceed as follows: 
dT iv 2fi{x) 



dx 

m 

dx 



B(x) 
1 



*(x) 

2/i(xQ 
B(x') 



^(x')dx' 



(20) 
(21) 



If a > then B{x) has an interior zero at x( a ) an d we integrate separately 
on (0,x(a)) and (x(a),2vr): 

*(cc) Jo B(x') I aX > 

_J_ « 2/^x0 ^/ ,x , / 



<9x 



(s) 



x < x(a) 
x > x(a) 



(22) 



Finally Tj(0) is given by a second integration: 



Ti(0) 



2 IT 



dx 



{x) dx 



(23) 



Here, the argument of the exponential function in ^f(x) can be any 



antiderivative of 



2A(x) 



B(x)-_ 



^(0) (and ^(x(a)), if a > 0) is zero because 
the adjoint equation (|14|) has an irregular singularity at x = (and at 

0+ 



x = x(«)); consequently |^ 



oo and j x dx 2A ^ 



B(x) 



-oo as x 



(and as x — > x( a ) )• This accounts for the difference between (23) and, for 
example, Eqn. (5.2.157) in [20] [] Note that, also because Z(0) = and 
<jj > 0, x = is an entrance boundary, so any exit must take place at 2tt 
(i.e. an exit from the interval (0,2-7r) is equivalent to a spike). 

For the class of PRCs that we consider, the antiderivative in (19) can be 
evaluated symbolically, therefore ^f(x) can be evaluated analytically. The 
integrals in (21 23) must be evaluated by numerical quadrature. 



We next discuss the integrability of (21 ) and (|22j). First, we consider the 
case of the theta neuron (a = 0), for which B{x) has only two zeros; at and 
2-7T. We can confirm the integrability by checking the following conditions: 

^'2/(x'), 



lim 

x^0+ 



lim 

X— *2n~ 



B(x r 



-ty{x)dx 
lim fy(x) 

x^0+ 



2f(x>) 



B(x' 



lim 

x^2w 



fy(x')dx 
#(x) 






±oo 
oo 



(24) 
(25) 

(26) 
(27) 



1 The expression we obtain is identical to the evaluation of this expression in the case 
of a reflecting boundary at 0. 
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If these are satisfied then by l'Hospital's rule, 



m 

dx 



1 

fy(x) 



c 2fi(x') 
B(x') 



^(x')dx' 



is finite at the endpoints; in fact 



hm — — 

x^O OX 

hm — — 

x^2n OX 



fi{x) 



lim . 
x->o A(x) 



fi{x) 



lim , 

x^2tt A(x) 



(28) 

(29) 
(30) 



and we can use a quadrature method that can handle integrable singularities. 
For a > 0, B{x) has one interior zero, dividing the domain into two intervals 



on which (14) has irregular singularities at each end. ^ must be computed 
separately on each interval. Again we find that ^* is finite on each interval, 



permitting computation of (22) with a standard quadrature routine. 



Finally, we compute the derivative of the firing rate with respect to a 



DC input; i.e. ^ for the system 



d9 = Ludt + Z(6)(fj,dt + aodW t ), G [0,2vr) 



(31) 



which is equivalent to the Ito SDE 



d9 



a 



{lo + y Z(0)Z'(6) + fiZ(9)) dt + aZ(9)dW u 



#£[0,2tt). (32) 



We wish to differentiate v with respect to \i and evaluate at ji = 0. According 
to equation (13) we must evaluate ^-(0,yu) where the drift term A(x,fi) = 



u) + \Z{x)Z \x) + fxZ(x) is a function of both x and fi. In general, T\ 
T2 and are also functions of two arguments (e.g. Ti(x,/x)) and we have 
indicated the arguments where needed for clarity. The notation 
to differentiation with respect to the first argument. 
We first consider the case a 



will refer 



0. Rewriting (23), we have 



1 



27T *(<E, AO B ( X 



—^(x, fj,)dx dx 



where 



^f(x, n) = exp 



2^ ; 

B(x') 



(33) 
(34) 
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Therefore 
dTi 



dp 



dp 



-2 tf(x» 



dx'dx 



nJo B{X>) #(x,/i) 



'2tt JO 

We use the relationship 



S(x') 



*(x,/i) 2 
2 &4(y,//) 

o ^(y) ^ 



dy 



dx'd$b) 



(36) 



to find that 
dTi 



dp 



(o,p) 



rx 



V(x',p) 



B{x>)V{x,p) J x , B{y) dp 



2 dA(y,p) 



dy dx' dx 



r-x r-y 



2 2 ftAfo./O 



„./„ ;„ BM »(„) B (,) a„ (37) 

assuming that the order of integration over x' and y can be switched. By 
Tonelli's theorem, this is valid if the integrand is single-signed. For a = 0, 
the integrand is always nonnegative (note that §4(#) = Z(x) and that B(x) 
is nonnegative). 



Next, we establish the integrability of the expression (37). Rewriting, 
we have 



dl\ 
dp 



(0,p) 



1 



2 OA 
2tt^{x,h)J B(y) dp 
1 



*(y,At) 



As j£(x) = Z(x) and are bounded, we can use the same conditions for 
integrability unchanged from d24 
A parallel derivation to (|33|- 
we have 



37| can be made when a > 0. In this case 



dTi 
dp 



(o,p) 



(38) 



where we have used the boundary condition Ti(2ir, p) = 0. The integrand 
is given by 



d dTi 
dp dx 



(x,p) 



x rv 2 y(x' :l i) 2 dA ^'^ (ix'dy 



Jo Jo B{x') *(x,At) B(y) dn 

rx ry 2 ttfoV) 2 dA(y,fi) 7 / » ^ / \ 



X < x(") 



(39) 
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Here, the switch in the order of integration is justified by the fact that the 
integrand ^ is positive on (0,x(a)) and negative on (x(a),27r), while the 
range of integration in Eqn. ( |39[ ) is always restricted to lie in one region or 
the other. 

3.3 Patterns of correlation transfer over long timescales 




Figure 2: Firing rate (left) and CV (right) of the theta model neuron over 
a range of parameters uo and a. Spike trains are illustrated for four sets 
of (to, a) values (see white squares), notably mean-driven (high u, low a - 
bottom spike train) and fluctuation-driven (low u>, high a - top spike train). 
Level sets of a are plotted (see text). 



Having derived and shown how to evaluate formula (10), we next use it 
to evaluate spike count correlations p. The results are formally valid in the 
limits of timescale T —* oo and input correlation c — ► 0, so we also conduct 
Monte Carlo simulations to reveal behavior of px for large but finite T and 
intermediate input correlation up to c = 0.3, and to test the applicability of 
our formula in these regimes. 

We compute these spike correlations over a range of u and a values that 
explores a full dynamical regime of the model. By this we mean that the 
values we use span from dominantly mean-driven firing ( e.g. co = 2.5, 
a = 0.4 at lower- right white square in Fig. [2|, to dominantly fluctuation- 
driven firing (e.g. oj = 0.4, a = 2.4 at upper-left white square), and all 
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Figure 3: (Left) Susceptibility for a = (top), a = 0.5 (middle), and 
a = 1 (bottom) over a range of parameters u and cr; lines are level sets of 
a = o/yfuj. (Right) px vs. c from Monte Carlo simulations for a = (top), 
a = 0.5, and a = 1 (bottom). Specific (w,cr) values shown as in Figure [2} 
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intermediate possibilities. 

The resulting output firing rate v, computed via Eqn. ( |11[ ), ranges from 
to 0.9 (measured in spikes per time unit); see Fig. |] (left panel). Note 
that v increases strongly with co and only weakly with a (see Sec. 3.4). The 
CV (via Eqn. (|12[)) ranges from to 0.55 (Fig. [2] right panel). We were 
unable to reach higher values of CV even with a much expanded range of 
a. By using a time change in the equations, CV can be seen to depend on 
the input parameters uj and a only through the relationship a = a / ' y/uj (see 
later in this section) and is therefore invariant on level curves of this ratio. 
As Fig. [2] shows, CV increases with this ratio. 

We first fix a moderately long time window T = 32, corresponding to 
1.6—16 interspike intervals for the parameter range at hand, and compute px 
from Monte Carlo simulations for a range of correlation strengths c £ [0, 0.3]; 
see Fig. [3] (right column). We see that, for Type I models, px is close to 
linear in this range of c, as for the linear response theory. For Type II 
models, linearity holds over a decreased range of c. Additionally, note that 
the limiting formula Eqn. (10) for p gives a close approximation to pt=32 
for Type I models in more fluctuation-driven regimes. The approximation 
is worse for dominantly mean-driven firing, and for Type II models, but the 
trend that correlations are lower for Type II models is correctly predicted 



by Eqn. 10 



Next, we discuss the trends in p predicted by the linear response theory. 
Two findings stand out in the left hand panels of Fig. |3j First, values of 
S(u>,a) (and hence p S c) are much larger for Type I (a = 0) than for 
Type II (a = 1) models. Second, S(uj,a) is nearly constant as the input 
mean and standard deviation u and a vary over a wide range, for both the 
Type I and Type II models. In sum, Type I models transfer « 66% of 
their input correlations into spike correlations over long timescales; Type II 
models transfer none of these input correlations, producing long-timescale 
spike counts that are uncorrelated. Additionally, an intermediate model 
(a = 1/2) transfers an intermediate level of correlations, and these levels do 
depend on ui and a. We will provide a partial explanation for overall trends 
in correlation transfer with a in Section 13.41 

Fig. [4] provides an alternative view of these results, by plotting the corre- 
lation gain S vs. the firing rate and CV that are evoked by input parameters 
drawn from the whole range of uj, a. First, note that S does not vary with 
firing rate for the Type I or Type II models, as expected from the previous 
plots. For the intermediate (a = 1/2) model, S does not display a clear 
functional relationship with firing rate, but there is such a relationship with 
CV (Fig. |4^b)). We note that all of these findings for the phase oscilla- 
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Figure 4: Susceptibility vs. v (a) and susceptibility vs. CV (b) for a = 
(light gray), 0.5 (medium gray), 1 (black). 



tors under study are in contrast to the behavior of linear integrate and fire 
neurons, which produce a strongly increasing, nearly functional relationship 
with firing rate |14|, |4"H] ; we revisit this point in the discussion. 

Now, we discuss a scaling relationship for the underlying equations that 
simplifies the parameter dependence and helps to explain the plots of S vs. 
v and S vs. CV. This symmetry (which was noted in [30j for the quadratic 
integrate and fire model), allows us to reduce the free parameters u,a to 
one parameter a = a/y/uj. The stochastic differential equation 

dd t = (u+°^Z{0)Z'{0) + iiZ{6)\dt + <TZ(6)dWt 
becomes, under a time change r = ut, 



d9 T = (l + —Z{d)Z'(0) + -Z(9)] dr + ^=Z(6)dW T 

\ 2i0 LO J \JUJ 



(40) 



~ 2 \ 

1 + —Z(e)Z'(0) + -Z{6) ) dr + aZ(6)dW T (41) 
2 oj J 

Each exit time must scale identically under this transformation, so that the 
exit time moments scale as 

v 



i 



uTi -» v = - (42) 



LO 



T 2 = to 2 T 2 (43) 
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and therefore the CV = yT^ — T^/Ti is invariant under the time change. 
Thus, the CV, which is computed with fi = 0, depends only on a, not on uj 
and a separately. Now consider the two sets of parameters (cj, a) and (1, cr), 
such that a = a/y/uj. According to (41), the effect of \i is scaled by 1/uj so 
that 



dv 
d/i 



Therefore, the susceptibility 



S(u, a) 



1 dv 

lo d/i 

°\%? 
CV 2 v 

cv 2 * 

o I dv \ 2 

w CV 2 v 
5(1, Si) 



(44) 

(45) 
(46) 

(47) 
(48) 



is invariant under the change of parameters (u>, a) — > (1, a/ y/u>)) exactly the 
same change of parameters under which the CV is conserved. Therefore, 
if there is a single contour for each value of CV (i.e., there are no discon- 
nected contours), we expect that S will be a function of CV, as in Fig.|4^b). 
Conversely, note that the firing rate will vary widely over any particular 
a contour, so that many different firing rates can be expected to yield the 
same S; unless S is constant, we expect a given firing rate to be associated 
with a range of S values, as also seen in Fig. [4] (a). 

Finally, in Fig. [5j we demonstrate that the behavior of S seems to 
roughly "interpolate" between the a = 0, 1/2, 1 cases considered here as 
it varies over it whole range of a. S typically decreases as the total noise 



variance a increases. 



3.4 Analytical arguments for Type I vs. 
in correlation gain 5* 



Type II difference 



We now seek to explain the drop in S(uj,a) (Eqn. (10)) as a increases, 
ranging from Type I (a = 0) to Type II (a = 1) PRCs. There are two 
tractable limits in which explicit calculations can be performed. First, we 
show that S(oj, a) = for "purely" Type II models, for any values of u) and 
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S(1,a) 




Figure 5: 5(1, a) for a continuum of models a £ [0, 1]. 



a. Next, for arbitrary a, we derive an S(u), a) valid to first order in a 2 - this 
reveals a monotonic decrease in 5 with a and gives a good description of 
correlation transfer even for a ~ 1. Finally, we given an intuitive argument 
to buttress these calculations. 



5 = for purely Type II models 

We start with a = 1. It is straightforward to see that dv/dfi = for any 
u, a. Recall that dv/d^i is given by the integral of a function over the 



interval (0,2-7r), in Eqn. (38); we will show this function integrates to zero. 

7T, 



Rewriting Eqn. (39), using x(l 



d dTi 
dfi dx 



4/;s®|(w) $ (w) x i(w)*. x 



< TT 
> TT. 



(49) 



fy(x, n),B(x), and ^(x,^) are each 7r-periodic (i.e. f(x) = f(x + ir)). 
Z{x) = — sin(x), however, is anti-periodic (/(x) = — fix + tt)). Then 



M 
dfi 



d dT x 
dfj, dx 



d dTi 
d\x dx 



(x + vr,/i) 



(50) 



and integrating y^ki x ^)j over a period x G (0, 27r) yields zero. 
Plugging this into Eqn. ( |10[ ), we see that 5(a>, a) = for models with Type 
II PRCs ZiO) = — sin(#) (z/ and CV are always nonzero in this paper). 
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Evaluation of S(u, a) in the limit a — ► 

We next derive an analytical expression for S(u,a) in the limit a —* 0. 
We will show that each relevant term Ti(0), 72(0), and ^^-(0) admits an 
asymptotic expansion in the small parameter a 2 . We compute these terms 
explicitly and combine them to get the first term of the associated expansion 
for S(uj, a). 

If we examine the integral in (21 ) and the inner integral of (35 ) for a = 0, 



or the integrals in (22), and (39) for a > 0, we see that we can write each 



of them in the form 
1 



z(y) 



f(x') exp ( —y(F(x') — F(y)) ) dx' 



(51) 



where F(x') is strictly increasing on (a, y). a may be either or x Q , depend- 
ing on the circumstance. This integral admits an asymptotic expansion in 
a 2 , with successive terms essentially given by integrating by parts and retain- 
ing only the contribution from the end of the interval where F{x')—F{y) = 0. 
In the case of T\(x) we have 



1 



In the case of Tz(x) we have 



and for dvjd\x 



Z(y) 

Ti(y) 
z{y) 

dT x 

dx 



(»)■ 



In each case F(x') is given by the antiderivative of z ^,y2 because this 
is a positive function on (0, x( a )) an( ^ (x( a )^ 7r )i F(x') must clearly be 
increasing on these intervals. The integral 



f(x') exp (A0(x')) dx' 



(52) 



admits the following expansion for A — > 00, if (/>' > on [a, b] and / meets 
certain conditions (see, for example, [7]): 



/(A) 



v (-1)" 

n=0 



exp (X(p(x)) 



1 d 

4>'(x) dx 



n fix) 



(j)'{x) 



(53) 
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If <f)(x) — ► — oo as x — > a, as is the case in our use of the formula, only the 
right-hand endpoint makes a contribution to the integral. 

Substituting this contribution into the outer integral and evaluating 
these quantities to the required order, we find 

T l(0 ) = — + 0(a i ) (54) 

T a (0) = 4 ^2+° 2 (^ } (3-6a + 4a 2 )) + 0(a 4 ) (55) 



^n(O) = --(l-«) + 0(a 4 ) (56) 
a/i a; 

In passing, we note that the firing rate gain, dv/dfj,, is given by 

dv 1 dT\ , . , 

5 - w* 1 " (57) 

= ^ + 0(<r 4 ) (58) 
Putting these results together we see that 

s <^> = »-fa°L +o ' ga » <59) 

It can be readily checked that this function decreases monotonically from a 
value of 2/3 at a = 0, to a value of at a = 1. While this calculation is in 
the limit <r — » + , it in fact remains a good approximation for moderate a, in 
fact, even for a ~ to. Figure |6] shows the limiting value as well as computed 
S values at small (relative to uj = 1; a = 0.2) and moderate (<r = 1) values of 
a. The limiting value remains a good approximation throughout this range. 

An argument for general PRCs 

We close this by section by noting that, for arbitrary PRCs z(9) and small 

0", 

(■2-K 



dv/dn oc / Z(8)d6 . 



o 

The calculations showing this are in the appendix. While S(oj, a) = - ^^Jfjf 1 
has other terms that depend on the PRC, this calculation does suggest that 
S is likely to be smaller for PRCs with lower means in general, and lends 
some intuition to what drives the decrease in S with a. 
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a 

Figure 6: lim a ^ + S(l, a) for a continuum of models a 6 [0,1] (black 
dashed). This is a good approximation for small (a = 0.2; dark gray) and 
moderate (a = 1; light gray) values of a. 

4 Correlation over shorter timescales 

Figure [7] shows px computed for a sequence of finite time windows T. We 
can characterize the non-dimensional T in terms of its length in terms of a 
typical interspike interval (ISI) of the oscillator. The time window T = 1 
varies from 0.05 - 0.5 ISI, roughly, from left to right; the time window T = 32 
varies from 1.6 — 16 ISI. 

We see a striking dependence of transferred correlations on T. p%2 is 
larger for Type I than for Type II oscillators for most parameter values, 
consistent with our long time results. 

However, p\{oj,a) is smaller for Type I than for Type II for 95 % of 
parameter pairs (u>,cr). This is consistent with recent results on the response 
of phase oscillators to correlated noise [311 [19]; in particular, [31 j study the 
distribution of phase difference A9 = 81 — 82 between two oscillators driven 
by common noise and find that the probability that A8 = is greater for 
Type II than for Type I. As noted in the Introduction, this metric can be 
shown to have a direct relationship with our px as T — > 0. To summarize, 
the "switch" in St(u,ct) (from higher correlation gain in Type II models 
to higher correlation gain in Type I) occurs at time scales T over which 
each cell fires several spikes; such timescales are biologically relevant, as we 
discuss in the next section. 
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T = 1 T = 2 T = 8 T = 32 




Figure 7: pt{u, <t) as measured from Monte Carlo simulations for an increas- 
ing sequence of T, for simulations with a = (top row), a = 0.5 (middle) 
and a = 1 (bottom). The fraction of common variance is c = 0.1; therefore, 
to recover the approximate Sr(w, o) ~ pt(u},<j)/c, multiply by 10. The u; 
and a axes of each plot are the same. As T increases, St approaches the 
T — > oo limit illustrated in Fig. [3} 

5 Summary and discussion 

We asked how correlated input currents are transferred into correlated spike 
trains in a class of nonlinear phase models that are generic reductions of 
neural oscillators. Linear response methods, asymptotics, and Monte Carlo 
simulations gave the following answers: 

1. Over long timescales, Type I oscillators transfer 66% of incoming cur- 
rent correlations into correlated spike counts, while Type II oscillators 
transfer almost none of their input correlations into spike count cor- 
relations. Models with intermediate phase response curves transfer 
intermediate levels of correlations. 

2. Over long timescales, correlation transfer in Type I and Type II models 
is independent of the rate and coefficient of variation (CV) of spiking. 
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For intermediate models, correlation transfer decreases with CV and 
shows no clear dependence on rate. 

3. That there is a timescale T beneath which these results reverse: Type 
II neurons become more efficient at transferring correlations than Type 
I, there is an increasing dependence of correlation transfer on spike 
rate, and the strong dependence on CV weakens. 

We note that results (1) and (2) are highly distinct from findings for the leaky 
integrate-and-fire neuron model, for which up to 90-100% of correlations 
are transferred over long timescales, with this level depending strongly on 
firing rate but very weakly on CV |14[ |4"8] . This demonstrates a strong 
role for subthreshold nonlinearities in determining correlation transfer in 
the oscillatory regime (as seen for the quadratic integrate-and-fire model 
in 08]). 

What timescales of spike count correlation actually matter in a given 
application? This depends on the circuit that is "downstream" of the pair 
(or, similarly, layer) of neural oscillators that we have studied in this pa- 
per; in other words, on what system is reading out the neurons we study 
here. Clearly, different neurons and networks are sensitive to input fluc- 
tuations over widely varying timescales. For example, some single neurons 
and circuits can respond only to events in which many of the cells that pro- 
vide inputs spike nearly simultaneously. This is the coincidence detector 
mode of operation (cf. jlQ] and references therin), and can result from fast 
membrane time constants (as occur in high-conductance states [HI]); circuit 
mechanisms, such as feed-forward inhibition [37], can also play a role. For 
such systems, short-timescale correlations among upstream cells are relevant 

- small window lengths T. On the opposite extreme, networks operating as 
neural integrators will accumulate inputs over arbitrarily long timescales 
(see [10] and references therein); in this case, spike-time correlations over 
large windows T are reflected in circuit activity. In general there is a range 
of possible behaviors, and the timescales over which inputs are integrated 
can differ among various components of a network, among components of an 
individual cell [S7\ I40j. or among different times in a cell lifetime, depending 
on background input characteristics |16j . 

One domain in which different levels of correlation transfer - and dif- 
ferent dependences of this correlation transfer on neurons' operating ranges 

- can have a strong effect is the population coding of sensory stimuli. For 
example, if neurons are read out over long timescales, then Type I vs Type 
II populations offer a choice between relatively high and low levels of corre- 
lation across the population. Depending on heterogeneity of the population 
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response to the stimulus at hand, one or the other of these choices can 
yield dramatically greater (Fisher) information about the encoded stimu- 
lus [5U CD H2] • The opposite choice of neuron type would be preferred for 
readout over short timescales, where trends in correlation transfer reverse. 
Beyond averaged levels of correlation, a separate question that can affect 
encoding is whether correlations depend on the stimulus |344 [23]. A nat- 
ural way that this can occur is when correlations depend on the evoked 
rate or CV of firing. We demonstrate that such dependencies are present in 
phase models over short timescales and, over long timescales, that they are 
present in intermediate but not "purely" Type-I or Type-II models. Once 
again, depending on details of stimulus encoding, these dependencies can 
either enhance or degrade encoding. Overall, the picture that emerges is 
that correlation transfer is another factor to consider in asking which non- 
linearities allow neuron models to best encode stimuli, and that there will 
be different answers for different stimuli. 

In closing, we note that we have studied only simple (but widely-used) 
one-dimensional neural models here. However, preliminary simulations sug- 
gest that the trends for correlation transfer found here also hold in some 
standard Type-I vs. Type-II conductance-based neuron models [38 . The 
situation is more complex, as global features of the neural dynamics can be 
involved, and will be explored in future work. 
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6 Appendix 



In this appendix we give some more details about calculation of 3i(0), ^2(0), 



and Trr(O) for specific values of a 



6.1 a = 1 and numerical details 

For Z(x) = — sin(x) (a = 1), 



^2 

A(x) = uj + — sin(x) cos(x) (60) 

B{x) = a 2 [sin(x)] 2 (61) 

are periodic on [0, tt] with B{x) = at 0(27r) and tt. 

fy(x) is defined as an anti-derivative as in Eqn. (19); we can compute 
this symbolically yielding 

2lo 

fy(x) = sin(x) exp( ~ cot(x)). (62) 

a A 

We can check that \\m x ^Q+ ^(x) = 0, linx,,^,,— ^f(x) = 00, and that ty/B is 
not integrable at either or tt. 

For the interval tt to 2tt we use the fact that A{x) and B{x) are ^-periodic; 
repeats on this second interval; that is 

^f(x) = I sin(x)| exp f ^ cot (a;) J (63) 

is an expression that is valid everywhere Z{6) 7^ 0. 

To compute values of on a uniform mesh, we use an adaptive quadra- 
ture routine to evaluate Eqns. (22); in particular Simpson's rule imple- 
mented via the MATLAB routine quad. We have already demonstrated 
the integrability of Eqns. (22) in £3.2 Ti(0) is then computed using Simp- 
son's 3-point rule. 

To compute values of we use adaptive quadrature as well, with the 
caveat that -S^ is evaluated in between mesh points by linear interpolation. 

6.2 a = 

For the theta model, Z{x) = 1 — cos(x) or a = 0, 

A(x) = uj + — (1 — cos(x)) sin(x) (64) 
B{x) = a 2 [l - cos(x)] 2 (65) 
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are periodic on [0, 2?r] with B{x) = at 0(27r). Again we can integrate 
2A{x) / B{x) symbolically, finding 

h / m ( 2w (2 -cos(x))sin(x)\ lnnS 
* = |1 - cos(s)| exp (- ^ (1 _L(V J (66) 



We can verify Eqn. (24) as before. 



6.2.1 A small a argument for general PRC 

Here we replicate the small a value of df/dfj,, using a perturbation expansion 
valid for an arbitrary PRC. 

We consider the stationary densities p{6) and p{9) of two different pro- 
cesses, 

d6 = a[9)dt + b{9)dW t (67) 
d9 = a(e)dt + b(0)dW t (68) 



where 



2 

a{6) = uj+ <J ^Z{e)Z'{e) (69) 

2 

a(o) = u + ^-z(e)z'(e) + fiZ(e) (70) 

p{9) satisfies the stationary Fokker-Planck equation 

~ We {ap) + W* {hp) = (71) 

This can be integrated and the constant of integration is equal to the firing 
rate: 

d 

ap--^{bp) = v (72) 



and therefore 



2tt 
1 

2^r 







2k 

a{e)p{e)de (73) 
a{e)p{e)de (74) 
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The DC input response, dv/dfj,, can be computed by differentiating v with 
respect to fi and evaluating at \i = 0. Therefore we expand v as a series in 
(j, and take the first-order term. 



P 
v 



P + HPI + 0(fi 2 ) 
v + nvx + 0(/i 2 ) 



(75) 
(76) 



We find that 
dv 



1 

2^ 



2tt 



2tt 



p{9)Z{e)d9+ / pi(9)a(9)d6 



(77) 



We consider (77) more carefully. We will see that all but one term are 
multiplied by a 2 ; for small <r, the term that remains significant is the mean 
of the phase-resetting curve. 



We rewrite Eqn. (77) as follows: 
I r /"27T 

1 

2^ 



o 



(9) P i(9)d9 



2tt 



^+P(0)) Z{9)d9 + a 2 J^ l -Z{9)Z' 

l -Z{9)d9+ p{9)Z{9)d9 + a 2 [*" \z(9)Z'(9) Pl (9iaA 
!7r Jo Jo 2 



using the fact that pi(#) must average to so that p remains a probability 
density. We have also written p(9), the stationary density for the process 
with fx = 0, as the sum of a uniform density and a deviation p{9). 
Consider the process 



d9 



to + ~Z(9)Z'{9) ) dt + aZ(9)dW t 



The stationary density p and firing rate J satisfy 



co + ^Z(9)Z'(9))p(9) 



d I a 
89 



-Z 2 (9)p{9) 



J 



(80) 



(81) 



If a = then the stationary density is po = ^- and Jo = We expand p 
and J in powers of cr 2 ; 



p(#) = Po+P 

= p + a 2 p 1 + 0(<7 4 ) 
J = J + <7 2 Ji + 0(a 4 ) 



(82) 
(83) 
(84) 
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At O(l) we have 



as already stated. At 0(a) 



Po = ~ (85) 



= - ( Jx + —ZZ'dd] (87) 
u> \ 4tt J 

J\ is determined so that p is a probability density at any order: 

2-7T i /*2-7r 

j 1 = 2^Ji = -— / ZZ'd9 = 



o 47r Jo 

as ZZ' is the perfect derivative of a periodic function. For general order, we 
have 

Pn = l( J n-\zZ'Pn-l+^{^Pn-ty (89) 
Jn = ^ ZZ'^dO (90) 

where we no longer expect J n to be zero. 



Let's return to (79). The first integral is the mean of Z(6). The second 
integral, in our expansion, first appears at fourth-order in a. To see this, we 
examine 

f piZ(6)d9 = — [ Z 2 (6)Z'(0)dO (91) 
Jo 4vr Jo 

= (92) 



The third term is second-order in a. Thus the dominant term in ( |79| ) is the 
first one, and we have shown the desired result. 
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